ASA Dynamic Risk Predictions Workshop notes

Part 1: Introduction

  • Explicit outcomes: multiple longitudinal responses. time-to-events

  • Implicit outcomes: missing data, random visit times

  • Research question with AIDS dataset: how strong is the association between CD4 cell count and the risk of death?

  • Research question with PBC dataset: how strong is the association between bilirubin and the risk of death, can it be used for predictions of survival probabilities

When to use separate analysis (only concerns either survival/longitudinal): * Does treatment affect survival? * Are average longitudinal evolutions different between males and females?

When to use joint analysis: * Complex effect estimation: how strong is the association between the longitudinal evolution of CD4 cell counts and the hazard of death? * Time-varying covariates vs. endogenous * Handling implicit outcomes: focus on longitudinal outcomes with dropouts and random visit times

Part 2: The joint-model

Framework: * Use model to describe the evolution of the covariate * Fit longitudinal data into cox model diagramofjointmodelframework

  • Think about fitting the green line into the hazard function

  • Random effects explains all interdependencies in the data

    • Longitudinal outcome is independent of TTE outcome
    • Repeated measurements in longitudinal are independent of each other
  • For cox model in joint, requires full specification of joint distribution

    • Use spline as a flexible model for baseline hazard
  • jm() function assumes that the ordering of the patients is the same, and that the time unit is the same

  • Great explanation of function at 1:10

  • endogeneous: inside patients, such as bio markers

  • extrogeneous: outside patients, such as nurse caring for patient, hospital

Part 3: Functional forms

  • Instantaneous markers are not always optimal, sometimes other functional forms are needed to accurately model
    • Time dependent
    • Cumulative effects (can be weighted)
  • functional form is specified by functional_forms
    • area() calculates the cumulative effect
    • slope() calculates the time-dependent slope

Part 4: Dynamic predictions

Dynamic prediction: survival probabilities are dynamically updated as additional longitudinal information is recorded

Fitted model estimates conditional survival probabilities

Model accounts all measurements present and before, but not in the future. Obviously becomes more accurate with more measurements

Predictions at 2:05

Functional forms have an effect on predictions, so choose the optimal one according to a model selection problem

Model can predict using new data, does not get re-run

Accuracy for predictions: * Discrimination: classifying high-risk and low-risk patients * Interested in events occuring in a medically-relevant interval * sensitivity, AUC and ROC curve * Use function tvROC() * Calibration: how well the joint model can accurately predict future events * Use a calibration plot: calibration_plot() * Brier score: combining discrimination and calibration * Lower brier score is better * Calculate using tvBrier() * Take the square root of the brier score to find the average probability difference (less than 5% is good) * To measure predictive ability, validation sets are needed

Practical run-through

Dataset: prothro, Liver Cirrhosis dataset

library("JMbayes2")
## Loading required package: survival
## Loading required package: nlme
## Loading required package: GLMMadaptive
## Loading required package: splines
# linear mixed effects model
lmeFit <- lme(pro ~ time + time:treat, data = prothro, random = ~ time | id)

# cox model
CoxFit <- coxph(Surv(Time, death) ~ treat, data = prothros)

# joint-model
jointFit <- jm(CoxFit, lmeFit, time_var = "time") 
# add functional_forms argument to specify functional form

summary(jointFit)
## 
## Call:
## jm(Surv_object = CoxFit, Mixed_objects = lmeFit, time_var = "time")
## 
## Data Descriptives:
## Number of Groups: 488        Number of events: 292 (59.8%)
## Number of Observations:
##   pro: 2968
## 
##                  DIC     WAIC      LPML
## marginal    28075.50 28058.46 -14029.52
## conditional 34725.74 34439.72 -17563.87
## 
## Random-effects covariance matrix:
##                      
##        StdDev    Corr
## (Intr) 19.3032 (Intr)
## time   4.0524  0.0201
## 
## Survival Outcome:
##                    Mean  StDev    2.5%   97.5%      P   Rhat
## treatprednisone -0.1727 0.1614 -0.4970  0.1315 0.2838 1.0008
## value(pro)      -0.0418 0.0040 -0.0496 -0.0343 0.0000 1.0123
## 
## Longitudinal Outcome: pro (family = gaussian, link = identity)
##                         Mean  StDev    2.5%   97.5%      P   Rhat
## (Intercept)          73.3298 0.9817 71.3660 75.2442 0.0000 0.9999
## time                  0.4516 0.4738 -0.4656  1.3774 0.3356 1.0259
## time:treatprednisone -0.1292 0.6579 -1.4263  1.1372 0.8496 1.0099
## sigma                17.2312 0.2560 16.7351 17.7553 0.0000 1.0059
## 
## MCMC summary:
## chains: 3 
## iterations per chain: 3500 
## burn-in per chain: 500 
## thinning: 1 
## time: 8 sec
  • Rule of thumb: Rhat needs to be below 1.1. If not, increase the number of iterations for MCMC

Extract the data of patient 155

dataP155 <- prothro[prothro$id == 155, ]
dataP155$Time <- dataP155$death <- NULL # delete variable as we want to predict this

Note: if you want another longitudinal biomarker, need to make another linear model for it

Prediction plot

sfit <- predict(jointFit, newdata = dataP155[1, ], process = "event", 
                times = seq(0, 10, length.out = 51),
                return_newdata = TRUE)

sfit
##         id pro    time      treat start      stop event death  Time  pred_CIF
## 809    155  38 1.0e-06 prednisone     0 0.2710546     0     0 1e-06 0.0000000
## 809.1  155  38 2.0e-01 prednisone     0 0.2710546     0     0 1e-06 0.1134616
## 809.2  155  38 4.0e-01 prednisone     0 0.2710546     0     0 1e-06 0.1989687
## 809.3  155  38 6.0e-01 prednisone     0 0.2710546     0     0 1e-06 0.2663289
## 809.4  155  38 8.0e-01 prednisone     0 0.2710546     0     0 1e-06 0.3212240
## 809.5  155  38 1.0e+00 prednisone     0 0.2710546     0     0 1e-06 0.3671941
## 809.6  155  38 1.2e+00 prednisone     0 0.2710546     0     0 1e-06 0.4065749
## 809.7  155  38 1.4e+00 prednisone     0 0.2710546     0     0 1e-06 0.4409817
## 809.8  155  38 1.6e+00 prednisone     0 0.2710546     0     0 1e-06 0.4715728
## 809.9  155  38 1.8e+00 prednisone     0 0.2710546     0     0 1e-06 0.4991285
## 809.10 155  38 2.0e+00 prednisone     0 0.2710546     0     0 1e-06 0.5241815
## 809.11 155  38 2.2e+00 prednisone     0 0.2710546     0     0 1e-06 0.5471486
## 809.12 155  38 2.4e+00 prednisone     0 0.2710546     0     0 1e-06 0.5683628
## 809.13 155  38 2.6e+00 prednisone     0 0.2710546     0     0 1e-06 0.5880936
## 809.14 155  38 2.8e+00 prednisone     0 0.2710546     0     0 1e-06 0.6065617
## 809.15 155  38 3.0e+00 prednisone     0 0.2710546     0     0 1e-06 0.6239497
## 809.16 155  38 3.2e+00 prednisone     0 0.2710546     0     0 1e-06 0.6404002
## 809.17 155  38 3.4e+00 prednisone     0 0.2710546     0     0 1e-06 0.6560065
## 809.18 155  38 3.6e+00 prednisone     0 0.2710546     0     0 1e-06 0.6708431
## 809.19 155  38 3.8e+00 prednisone     0 0.2710546     0     0 1e-06 0.6849745
## 809.20 155  38 4.0e+00 prednisone     0 0.2710546     0     0 1e-06 0.6984572
## 809.21 155  38 4.2e+00 prednisone     0 0.2710546     0     0 1e-06 0.7113425
## 809.22 155  38 4.4e+00 prednisone     0 0.2710546     0     0 1e-06 0.7236772
## 809.23 155  38 4.6e+00 prednisone     0 0.2710546     0     0 1e-06 0.7355044
## 809.24 155  38 4.8e+00 prednisone     0 0.2710546     0     0 1e-06 0.7468619
## 809.25 155  38 5.0e+00 prednisone     0 0.2710546     0     0 1e-06 0.7577816
## 809.26 155  38 5.2e+00 prednisone     0 0.2710546     0     0 1e-06 0.7682922
## 809.27 155  38 5.4e+00 prednisone     0 0.2710546     0     0 1e-06 0.7784205
## 809.28 155  38 5.6e+00 prednisone     0 0.2710546     0     0 1e-06 0.7881921
## 809.29 155  38 5.8e+00 prednisone     0 0.2710546     0     0 1e-06 0.7976317
## 809.30 155  38 6.0e+00 prednisone     0 0.2710546     0     0 1e-06 0.8067629
## 809.31 155  38 6.2e+00 prednisone     0 0.2710546     0     0 1e-06 0.8155952
## 809.32 155  38 6.4e+00 prednisone     0 0.2710546     0     0 1e-06 0.8241113
## 809.33 155  38 6.6e+00 prednisone     0 0.2710546     0     0 1e-06 0.8322932
## 809.34 155  38 6.8e+00 prednisone     0 0.2710546     0     0 1e-06 0.8401270
## 809.35 155  38 7.0e+00 prednisone     0 0.2710546     0     0 1e-06 0.8476036
## 809.36 155  38 7.2e+00 prednisone     0 0.2710546     0     0 1e-06 0.8547191
## 809.37 155  38 7.4e+00 prednisone     0 0.2710546     0     0 1e-06 0.8614750
## 809.38 155  38 7.6e+00 prednisone     0 0.2710546     0     0 1e-06 0.8678779
## 809.39 155  38 7.8e+00 prednisone     0 0.2710546     0     0 1e-06 0.8739400
## 809.40 155  38 8.0e+00 prednisone     0 0.2710546     0     0 1e-06 0.8796751
## 809.41 155  38 8.2e+00 prednisone     0 0.2710546     0     0 1e-06 0.8850968
## 809.42 155  38 8.4e+00 prednisone     0 0.2710546     0     0 1e-06 0.8902190
## 809.43 155  38 8.6e+00 prednisone     0 0.2710546     0     0 1e-06 0.8950558
## 809.44 155  38 8.8e+00 prednisone     0 0.2710546     0     0 1e-06 0.8996219
## 809.45 155  38 9.0e+00 prednisone     0 0.2710546     0     0 1e-06 0.9039324
## 809.46 155  38 9.2e+00 prednisone     0 0.2710546     0     0 1e-06 0.9080037
## 809.47 155  38 9.4e+00 prednisone     0 0.2710546     0     0 1e-06 0.9118521
## 809.48 155  38 9.6e+00 prednisone     0 0.2710546     0     0 1e-06 0.9154917
## 809.49 155  38 9.8e+00 prednisone     0 0.2710546     0     0 1e-06 0.9189352
## 809.50 155  38 1.0e+01 prednisone     0 0.2710546     0     0 1e-06 0.9221942
##           low_CIF   upp_CIF
## 809    0.00000000 0.0000000
## 809.1  0.03717772 0.3141227
## 809.2  0.06794302 0.5199672
## 809.3  0.09523238 0.6602401
## 809.4  0.12073251 0.7536370
## 809.5  0.14155980 0.8112603
## 809.6  0.15629531 0.8578796
## 809.7  0.16938166 0.8918004
## 809.8  0.18643586 0.9169940
## 809.9  0.20281096 0.9375112
## 809.10 0.21801270 0.9530290
## 809.11 0.23221165 0.9648173
## 809.12 0.24556498 0.9738068
## 809.13 0.25829630 0.9802099
## 809.14 0.27019033 0.9843905
## 809.15 0.27960896 0.9876075
## 809.16 0.28866421 0.9901030
## 809.17 0.29756012 0.9920583
## 809.18 0.31031221 0.9935954
## 809.19 0.32447311 0.9948176
## 809.20 0.33417186 0.9958026
## 809.21 0.34374247 0.9966013
## 809.22 0.35310894 0.9972520
## 809.23 0.36248599 0.9977843
## 809.24 0.37620089 0.9982213
## 809.25 0.39025750 0.9985811
## 809.26 0.40469455 0.9988776
## 809.27 0.41472089 0.9991215
## 809.28 0.42350592 0.9993579
## 809.29 0.43238794 0.9995927
## 809.30 0.43703992 0.9996821
## 809.31 0.44134054 0.9998156
## 809.32 0.44544157 0.9999178
## 809.33 0.44937362 0.9999560
## 809.34 0.45684734 0.9999824
## 809.35 0.46577067 0.9999953
## 809.36 0.47491189 0.9999989
## 809.37 0.48492628 0.9999998
## 809.38 0.49512198 1.0000000
## 809.39 0.50533315 1.0000000
## 809.40 0.51556616 1.0000000
## 809.41 0.52582714 1.0000000
## 809.42 0.53642914 1.0000000
## 809.43 0.54636979 1.0000000
## 809.44 0.55662572 1.0000000
## 809.45 0.56691491 1.0000000
## 809.46 0.57725642 1.0000000
## 809.47 0.58769190 1.0000000
## 809.48 0.59826448 1.0000000
## 809.49 0.60901724 1.0000000
## 809.50 0.61999342 1.0000000

Plotting prediction plot

plot(sfit)

Predicting the longitudinal outcome

Lpred <- predict(jointFit, newdata = dataP155[1, ],
                 times = seq(0, 10, length.out = 51),
                 return_newdata = TRUE)

plot(Lpred)

Combining the plots

plot(Lpred, sfit)

Dynamic predictions

n <- nrow(dataP155)
for (i in seq_len(n)) {
  Spred <- predict(jointFit, newdata = dataP155[1:i, ],
                   process = "event", times = seq(0, 10, length.out = 51),
                   return_newdata = TRUE)
  Lpred <- predict(jointFit, newdata = dataP155[1:i, ],
                   times = seq(0, 10, length.out = 51),
                   return_newdata = TRUE)
  plot(Lpred,Spred)
}

ROC at 3 years with 1 year window

roc <- tvROC(jointFit, newdata = prothro, Tstart = 3, Dt = 1)
plot(roc)

AUC at 3 years with 1 year window

tvAUC(roc)
## 
##  Time-dependent AUC for the Joint Model jointFit
## 
## Estimated AUC: 0.6819
## At time: 4
## Using information up to time: 3 (229 subjects still at risk)

Calibration plot at 3 years with 1 year window

calibration_plot(jointFit, newdata = prothro, Tstart = 3, Dt = 1)

Brier score at 3 years with 1 year window

tvBrier(jointFit, newdata = prothro, Tstart = 3, Dt = 1)
## 
## Prediction Error for the Joint Model 'jointFit'
## 
## Estimated Brier score: 0.1015
## At time: 4
## For the 229 subjects at risk at time 3
## Number of subjects with an event in [3, 4): 27
## Number of subjects with a censored time in [3, 4): 14
## Accounting for censoring using model-based weights

Concluding remarks

  • Pay attention to:
    • Model longitudinal flexibly using splines
    • Consider functional form for joint model

Tuning longitudinal data through splines

Spline overview: https://rpubs.com/alecri/review_longitudinal

Data: Study of 162 girls on body fat percentage pre-menarcheal and post-menarcheal

# loading packages
library(labelled)   # labeling data
library(rstatix)    # summary statistics
## 
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
## 
##     filter
library(ggpubr)     # convenient summary statistics and plots
## Loading required package: ggplot2
library(GGally)     # advanced plot
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
library(car)        # useful for anova/wald test
## Loading required package: carData
library(Epi)        # easy getting CI for model coef/pred
library(lme4)       # linear mixed-effects models
## Loading required package: Matrix
## 
## Attaching package: 'lme4'
## The following object is masked from 'package:Epi':
## 
##     factorize
## The following object is masked from 'package:GLMMadaptive':
## 
##     negative.binomial
## The following object is masked from 'package:nlme':
## 
##     lmList
library(lmerTest)   # test for linear mixed-effects models
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(emmeans)    # marginal means
## 
## Attaching package: 'emmeans'
## The following object is masked from 'package:GGally':
## 
##     pigs
library(multcomp)   # CI for linear combinations of model coef
## Loading required package: mvtnorm
## Loading required package: TH.data
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:rstatix':
## 
##     select
## The following object is masked from 'package:JMbayes2':
## 
##     area
## The following object is masked from 'package:GLMMadaptive':
## 
##     negative.binomial
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
library(geepack)    # generalized estimating equations
library(ggeffects)  # marginal effects, adjusted predictions
library(gt)         # nice tables
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.2     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ lubridate 1.9.2     ✔ tibble    3.2.1
## ✔ purrr     1.0.1     ✔ tidyr     1.3.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::collapse() masks nlme::collapse()
## ✖ tidyr::expand()   masks Matrix::expand()
## ✖ dplyr::filter()   masks rstatix::filter(), stats::filter()
## ✖ dplyr::lag()      masks stats::lag()
## ✖ tidyr::pack()     masks Matrix::pack()
## ✖ dplyr::recode()   masks car::recode()
## ✖ dplyr::select()   masks MASS::select(), rstatix::select()
## ✖ purrr::some()     masks car::some()
## ✖ tidyr::unpack()   masks Matrix::unpack()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
# load data
load(url("http://alecri.github.io/downloads/data/bodyfat.RData"))

head(bodyfat)
## # A tibble: 6 × 6
##      id  agec  agem    time   pbf occasion
##   <dbl> <dbl> <dbl>   <dbl> <dbl>    <dbl>
## 1     1  9.32  13.2 -3.87    7.94        1
## 2     1 10.3   13.2 -2.86   15.6         2
## 3     1 11.2   13.2 -1.95   13.5         3
## 4     1 12.2   13.2 -1      23.2         4
## 5     1 13.2   13.2  0.0500 10.5         5
## 6     1 14.2   13.2  1.05   20.5         6
# trajectories overtime
ggplot(bodyfat, aes(agec, pbf)) +
  geom_line(aes(group = id)) +
  geom_smooth(se = FALSE, size = 2) +
  geom_vline(xintercept = 12.86, linetype = "dashed") + 
  labs(x = "Age, years", y = "Percent Body Fat")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

# fitting model
bodyfat <- bodyfat %>%
  mutate(timepost = pmax(time, 0))
lin_lspl0 <- lmer(pbf ~ time + timepost + (time + timepost | id), data = bodyfat)
summary(lin_lspl0)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: pbf ~ time + timepost + (time + timepost | id)
##    Data: bodyfat
## 
## REML criterion at convergence: 6062.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7742 -0.5900 -0.0359  0.5946  3.3798 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr       
##  id       (Intercept) 45.942   6.778               
##           time         1.631   1.277     0.29      
##           timepost     2.750   1.658    -0.54 -0.83
##  Residual              9.473   3.078               
## Number of obs: 1049, groups:  id, 162
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  21.3614     0.5646 161.5587  37.837  < 2e-16 ***
## time          0.4171     0.1572 108.4615   2.654  0.00915 ** 
## timepost      2.0471     0.2280 132.6814   8.980 2.32e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) time  
## time      0.351       
## timepost -0.515 -0.872
# slopes pre-menarcheal and post-menarcheal
K <- rbind(
  "population mean pre-menarcheal slope" = c(0, 1, 0),
  "population mean post-menarcheal slope" = c(0, 1, 1)
)
tidy(glht(lin_lspl0, linfct = K), conf.int = TRUE)
## # A tibble: 2 × 7
##   contrast          null.value estimate conf.low conf.high statistic adj.p.value
##   <chr>                  <dbl>    <dbl>    <dbl>     <dbl>     <dbl>       <dbl>
## 1 population mean …          0    0.417   0.0674     0.767      2.65      0.0156
## 2 population mean …          0    2.46    2.20       2.73      20.7       0
sid <- c(4, 16)
expand.grid(
  time = seq(-6, 4.1, 0.05),
  id = sid
) %>%
  mutate(timepost = pmax(time, 0)) %>%
  bind_cols(
    indiv_pred = predict(lin_lspl0, newdata = .),
    marg_pred = predict(lin_lspl0, newdata = ., re.form = ~ 0)
  ) %>%
ggplot(aes(time, indiv_pred, group = id, col = factor(id))) +
  geom_line(aes(y = marg_pred), col = "black", lwd = 1.5) +
  geom_line(aes(linetype = "BLUP")) + #empirical best linear unbiased prediction
  geom_point(data = filter(bodyfat, id %in% sid),
             aes(y = pbf, shape = factor(id), col = factor(id))) +
  guides(shape = "none") +
  labs(x = "Time relative to menarche, years", y = "Percent Body Fat",
       col = "Id", linetype = "Prediction")

Source: http://users.stat.umn.edu/~helwig/notes/smooth-spline-notes.html

Data: prestige of Canadian occupations from 1971 along with average income

# load data
library(car)

data(Prestige)
head(Prestige)
##                     education income women prestige census type
## gov.administrators      13.11  12351 11.16     68.8   1113 prof
## general.managers        12.26  25879  4.02     69.1   1130 prof
## accountants             12.77   9271 15.70     63.4   1171 prof
## purchasing.officers     11.42   8865  9.11     56.8   1175 prof
## chemists                14.62   8403 11.68     73.5   2111 prof
## physicists              15.64  11030  5.13     77.6   2113 prof
# plot data
plot(Prestige$income, Prestige$prestige,
     xlab = "Income", ylab = "Prestige")

# fit model
library(npreg)

mod.ss <- with(Prestige, ss(income, prestige))
mod.ss
## 
## Call:
## ss(x = income, y = prestige)
## 
## Smoothing Parameter  spar = 0.4923283   lambda = 0.0002148891
## Equivalent Degrees of Freedom (Df) 3.467512
## Penalized Criterion (RSS) 12118.05
## Generalized Cross-Validation (GCV) 127.3134
summary(mod.ss)
## 
## Call:
## ss(x = income, y = prestige)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.991  -7.727  -2.453   7.586  33.006 
## 
## Approx. Signif. of Parametric Effects:
##             Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)    63.08      2.228  28.310 0.000e+00 ***
## x              61.01      8.035   7.593 1.807e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Approx. Signif. of Nonparametric Effects:
##               Df Sum Sq Mean Sq F value    Pr(>F)    
## s(x)       1.468   2125    1448   11.78 0.0001596 ***
## Residuals 98.532  12118     123                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Residual standard error: 11.09 on 98.53 degrees of freedom
## Multiple R-squared:  0.5949,    Adjusted R-squared:  0.5845
## F-statistic: 57.35 on 2.468 and 98.53 DF,  p-value: <2e-16
plot(mod.ss, xlab = "Income", ylab = "Prestige")
with(Prestige, points(income, prestige))

Review on survival models

Source: https://www.geeksforgeeks.org/survival-analysis-in-r/

Survival analysis is the prediction of events at a specified time.

There are two methods of survival analysis in R: * Kaplan-Meier * Cox proportional hazard model

Kaplan-Meier

Used for censored data, not based on underlying probability distribution

library(survival)
lung
##     inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
## 1      3  306      2  74   1       1       90       100     1175      NA
## 2      3  455      2  68   1       0       90        90     1225      15
## 3      3 1010      1  56   1       0       90        90       NA      15
## 4      5  210      2  57   1       1       90        60     1150      11
## 5      1  883      2  60   1       0      100        90       NA       0
## 6     12 1022      1  74   1       1       50        80      513       0
## 7      7  310      2  68   2       2       70        60      384      10
## 8     11  361      2  71   2       2       60        80      538       1
## 9      1  218      2  53   1       1       70        80      825      16
## 10     7  166      2  61   1       2       70        70      271      34
## 11     6  170      2  57   1       1       80        80     1025      27
## 12    16  654      2  68   2       2       70        70       NA      23
## 13    11  728      2  68   2       1       90        90       NA       5
## 14    21   71      2  60   1      NA       60        70     1225      32
## 15    12  567      2  57   1       1       80        70     2600      60
## 16     1  144      2  67   1       1       80        90       NA      15
## 17    22  613      2  70   1       1       90       100     1150      -5
## 18    16  707      2  63   1       2       50        70     1025      22
## 19     1   61      2  56   2       2       60        60      238      10
## 20    21   88      2  57   1       1       90        80     1175      NA
## 21    11  301      2  67   1       1       80        80     1025      17
## 22     6   81      2  49   2       0      100        70     1175      -8
## 23    11  624      2  50   1       1       70        80       NA      16
## 24    15  371      2  58   1       0       90       100      975      13
## 25    12  394      2  72   1       0       90        80       NA       0
## 26    12  520      2  70   2       1       90        80      825       6
## 27     4  574      2  60   1       0      100       100     1025     -13
## 28    13  118      2  70   1       3       60        70     1075      20
## 29    13  390      2  53   1       1       80        70      875      -7
## 30     1   12      2  74   1       2       70        50      305      20
## 31    12  473      2  69   2       1       90        90     1025      -1
## 32     1   26      2  73   1       2       60        70      388      20
## 33     7  533      2  48   1       2       60        80       NA     -11
## 34    16  107      2  60   2       2       50        60      925     -15
## 35    12   53      2  61   1       2       70       100     1075      10
## 36     1  122      2  62   2       2       50        50     1025      NA
## 37    22  814      2  65   1       2       70        60      513      28
## 38    15  965      1  66   2       1       70        90      875       4
## 39     1   93      2  74   1       2       50        40     1225      24
## 40     1  731      2  64   2       1       80       100     1175      15
## 41     5  460      2  70   1       1       80        60      975      10
## 42    11  153      2  73   2       2       60        70     1075      11
## 43    10  433      2  59   2       0       90        90      363      27
## 44    12  145      2  60   2       2       70        60       NA      NA
## 45     7  583      2  68   1       1       60        70     1025       7
## 46     7   95      2  76   2       2       60        60      625     -24
## 47     1  303      2  74   1       0       90        70      463      30
## 48     3  519      2  63   1       1       80        70     1025      10
## 49    13  643      2  74   1       0       90        90     1425       2
## 50    22  765      2  50   2       1       90       100     1175       4
## 51     3  735      2  72   2       1       90        90       NA       9
## 52    12  189      2  63   1       0       80        70       NA       0
## 53    21   53      2  68   1       0       90       100     1025       0
## 54     1  246      2  58   1       0      100        90     1175       7
## 55     6  689      2  59   1       1       90        80     1300      15
## 56     1   65      2  62   1       0       90        80      725      NA
## 57     5    5      2  65   2       0      100        80      338       5
## 58    22  132      2  57   1       2       70        60       NA      18
## 59     3  687      2  58   2       1       80        80     1225      10
## 60     1  345      2  64   2       1       90        80     1075      -3
## 61    22  444      2  75   2       2       70        70      438       8
## 62    12  223      2  48   1       1       90        80     1300      68
## 63    21  175      2  73   1       1       80       100     1025      NA
## 64    11   60      2  65   2       1       90        80     1025       0
## 65     3  163      2  69   1       1       80        60     1125       0
## 66     3   65      2  68   1       2       70        50      825       8
## 67    16  208      2  67   2       2       70        NA      538       2
## 68     5  821      1  64   2       0       90        70     1025       3
## 69    22  428      2  68   1       0      100        80     1039       0
## 70     6  230      2  67   1       1       80       100      488      23
## 71    13  840      1  63   1       0       90        90     1175      -1
## 72     3  305      2  48   2       1       80        90      538      29
## 73     5   11      2  74   1       2       70       100     1175       0
## 74     2  132      2  40   1       1       80        80       NA       3
## 75    21  226      2  53   2       1       90        80      825       3
## 76    12  426      2  71   2       1       90        90     1075      19
## 77     1  705      2  51   2       0      100        80     1300       0
## 78     6  363      2  56   2       1       80        70     1225      -2
## 79     3   11      2  81   1       0       90        NA      731      15
## 80     1  176      2  73   1       0       90        70      169      30
## 81     4  791      2  59   1       0      100        80      768       5
## 82    13   95      2  55   1       1       70        90     1500      15
## 83    11  196      1  42   1       1       80        80     1425       8
## 84    21  167      2  44   2       1       80        90      588      -1
## 85    16  806      1  44   1       1       80        80     1025       1
## 86     6  284      2  71   1       1       80        90     1100      14
## 87    22  641      2  62   2       1       80        80     1150       1
## 88    21  147      2  61   1       0      100        90     1175       4
## 89    13  740      1  44   2       1       90        80      588      39
## 90     1  163      2  72   1       2       70        70      910       2
## 91    11  655      2  63   1       0      100        90      975      -1
## 92    22  239      2  70   1       1       80       100       NA      23
## 93     5   88      2  66   1       1       90        80      875       8
## 94    10  245      2  57   2       1       80        60      280      14
## 95     1  588      1  69   2       0      100        90       NA      13
## 96    12   30      2  72   1       2       80        60      288       7
## 97     3  179      2  69   1       1       80        80       NA      25
## 98    12  310      2  71   1       1       90       100       NA       0
## 99    11  477      2  64   1       1       90       100      910       0
## 100    3  166      2  70   2       0       90        70       NA      10
## 101    1  559      1  58   2       0      100       100      710      15
## 102    6  450      2  69   2       1       80        90     1175       3
## 103   13  364      2  56   1       1       70        80       NA       4
## 104    6  107      2  63   1       1       90        70       NA       0
## 105   13  177      2  59   1       2       50        NA       NA      32
## 106   12  156      2  66   1       1       80        90      875      14
## 107   26  529      1  54   2       1       80       100      975      -3
## 108    1   11      2  67   1       1       90        90      925      NA
## 109   21  429      2  55   1       1      100        80      975       5
## 110    3  351      2  75   2       2       60        50      925      11
## 111   13   15      2  69   1       0       90        70      575      10
## 112    1  181      2  44   1       1       80        90     1175       5
## 113   10  283      2  80   1       1       80       100     1030       6
## 114    3  201      2  75   2       0       90       100       NA       1
## 115    6  524      2  54   2       1       80       100       NA      15
## 116    1   13      2  76   1       2       70        70      413      20
## 117    3  212      2  49   1       2       70        60      675      20
## 118    1  524      2  68   1       2       60        70     1300      30
## 119   16  288      2  66   1       2       70        60      613      24
## 120   15  363      2  80   1       1       80        90      346      11
## 121   22  442      2  75   1       0       90        90       NA       0
## 122   26  199      2  60   2       2       70        80      675      10
## 123    3  550      2  69   2       1       70        80      910       0
## 124   11   54      2  72   1       2       60        60      768      -3
## 125    1  558      2  70   1       0       90        90     1025      17
## 126   22  207      2  66   1       1       80        80      925      20
## 127    7   92      2  50   1       1       80        60     1075      13
## 128   12   60      2  64   1       1       80        90      993       0
## 129   16  551      1  77   2       2       80        60      750      28
## 130   12  543      1  48   2       0       90        60       NA       4
## 131    4  293      2  59   2       1       80        80      925      52
## 132   16  202      2  53   1       1       80        80       NA      20
## 133    6  353      2  47   1       0      100        90     1225       5
## 134   13  511      1  55   2       1       80        70       NA      49
## 135    1  267      2  67   1       0       90        70      313       6
## 136   22  511      1  74   2       2       60        40       96      37
## 137   12  371      2  58   2       1       80        70       NA       0
## 138   13  387      2  56   1       2       80        60     1075      NA
## 139    1  457      2  54   1       1       90        90      975      -5
## 140    5  337      2  56   1       0      100       100     1500      15
## 141   21  201      2  73   2       2       70        60     1225     -16
## 142    3  404      1  74   1       1       80        70      413      38
## 143   26  222      2  76   1       2       70        70     1500       8
## 144    1   62      2  65   2       1       80        90     1075       0
## 145   11  458      1  57   1       1       80       100      513      30
## 146   26  356      1  53   2       1       90        90       NA       2
## 147   16  353      2  71   1       0      100        80      775       2
## 148   16  163      2  54   1       1       90        80     1225      13
## 149   12   31      2  82   1       0      100        90      413      27
## 150   13  340      2  59   2       0      100        90       NA       0
## 151   13  229      2  70   1       1       70        60     1175      -2
## 152   22  444      1  60   1       0       90       100       NA       7
## 153    5  315      1  62   2       0       90        90       NA       0
## 154   16  182      2  53   2       1       80        60       NA       4
## 155   32  156      2  55   1       2       70        30     1025      10
## 156   NA  329      2  69   1       2       70        80      713      20
## 157   26  364      1  68   2       1       90        90       NA       7
## 158    4  291      2  62   1       2       70        60      475      27
## 159   12  179      2  63   1       1       80        70      538      -2
## 160    1  376      1  56   2       1       80        90      825      17
## 161   32  384      1  62   2       0       90        90      588       8
## 162   10  268      2  44   2       1       90       100     2450       2
## 163   11  292      1  69   1       2       60        70     2450      36
## 164    6  142      2  63   1       1       90        80      875       2
## 165    7  413      1  64   1       1       80        70      413      16
## 166   16  266      1  57   2       0       90        90     1075       3
## 167   11  194      2  60   2       1       80        60       NA      33
## 168   21  320      2  46   1       0      100       100      860       4
## 169    6  181      2  61   1       1       90        90      730       0
## 170   12  285      2  65   1       0      100        90     1025       0
## 171   13  301      1  61   1       1       90       100      825       2
## 172    2  348      2  58   2       0       90        80     1225      10
## 173    2  197      2  56   1       1       90        60      768      37
## 174   16  382      1  43   2       0      100        90      338       6
## 175    1  303      1  53   1       1       90        80     1225      12
## 176   13  296      1  59   2       1       80       100     1025       0
## 177    1  180      2  56   1       2       60        80     1225      -2
## 178   13  186      2  55   2       1       80        70       NA      NA
## 179    1  145      2  53   2       1       80        90      588      13
## 180    7  269      1  74   2       0      100       100      588       0
## 181   13  300      1  60   1       0      100       100      975       5
## 182    1  284      1  39   1       0      100        90     1225      -5
## 183   16  350      2  66   2       0       90       100     1025      NA
## 184   32  272      1  65   2       1       80        90       NA      -1
## 185   12  292      1  51   2       0       90        80     1225       0
## 186   12  332      1  45   2       0       90       100      975       5
## 187    2  285      2  72   2       2       70        90      463      20
## 188    3  259      1  58   1       0       90        80     1300       8
## 189   15  110      2  64   1       1       80        60     1025      12
## 190   22  286      2  53   1       0       90        90     1225       8
## 191   16  270      2  72   1       1       80        90      488      14
## 192   16   81      2  52   1       2       60        70     1075      NA
## 193   12  131      2  50   1       1       90        80      513      NA
## 194    1  225      1  64   1       1       90        80      825      33
## 195   22  269      2  71   1       1       90        90     1300      -2
## 196   12  225      1  70   1       0      100       100     1175       6
## 197   32  243      1  63   2       1       80        90      825       0
## 198   21  279      1  64   1       1       90        90       NA       4
## 199    1  276      1  52   2       0      100        80      975       0
## 200   32  135      2  60   1       1       90        70     1275       0
## 201   15   79      2  64   2       1       90        90      488      37
## 202   22   59      2  73   1       1       60        60     2200       5
## 203   32  240      1  63   2       0       90       100     1025       0
## 204    3  202      1  50   2       0      100       100      635       1
## 205   26  235      1  63   2       0      100        90      413       0
## 206   33  105      2  62   1       2       NA        70       NA      NA
## 207    5  224      1  55   2       0       80        90       NA      23
## 208   13  239      2  50   2       2       60        60     1025      -3
## 209   21  237      1  69   1       1       80        70       NA      NA
## 210   33  173      1  59   2       1       90        80       NA      10
## 211    1  252      1  60   2       0      100        90      488      -2
## 212    6  221      1  67   1       1       80        70      413      23
## 213   15  185      1  69   1       1       90        70     1075       0
## 214   11   92      1  64   2       2       70       100       NA      31
## 215   11   13      2  65   1       1       80        90       NA      10
## 216   11  222      1  65   1       1       90        70     1025      18
## 217   13  192      1  41   2       1       90        80       NA     -10
## 218   21  183      2  76   1       2       80        60      825       7
## 219   11  211      1  70   2       2       70        30      131       3
## 220    2  175      1  57   2       0       80        80      725      11
## 221   22  197      1  67   1       1       80        90     1500       2
## 222   11  203      1  71   2       1       80        90     1025       0
## 223    1  116      2  76   1       1       80        80       NA       0
## 224    1  188      1  77   1       1       80        60       NA       3
## 225   13  191      1  39   1       0       90        90     2350      -5
## 226   32  105      1  75   2       2       60        70     1025       5
## 227    6  174      1  66   1       1       90       100     1075       1
## 228   22  177      1  58   2       1       80        90     1060       0

Fitting the model

survival_function = survfit(Surv(lung$time, lung$status == 2) ~ 1)
survival_function
## Call: survfit(formula = Surv(lung$time, lung$status == 2) ~ 1)
## 
##        n events median 0.95LCL 0.95UCL
## [1,] 228    165    310     285     363

Plotting the function

plot(survival_function, main = "Kaplan-Meier", xlab = "Number of days", 
     ylab = "Prob of survival")

Cox proportional hazard model

Uses the hazard function, which considers independent variables in regression.

Also does not assume an underlying probability, but assumes that the hazards are constant over time.

More useful to measure instantaneous risk of deaths, often delivers better results because more volatile with data and features.

Tends to drop sharper as time increases.

cox_mod <- coxph(Surv(lung$time, lung$status == 2) ~., data = lung)

summary(cox_mod)
## Call:
## coxph(formula = Surv(lung$time, lung$status == 2) ~ ., data = lung)
## 
##   n= 167, number of events= 120 
##    (61 observations deleted due to missingness)
## 
##                 coef  exp(coef)   se(coef)      z Pr(>|z|)    
## inst      -3.037e-02  9.701e-01  1.312e-02 -2.315 0.020619 *  
## age        1.281e-02  1.013e+00  1.194e-02  1.073 0.283403    
## sex       -5.666e-01  5.674e-01  2.014e-01 -2.814 0.004890 ** 
## ph.ecog    9.074e-01  2.478e+00  2.386e-01  3.803 0.000143 ***
## ph.karno   2.658e-02  1.027e+00  1.163e-02  2.286 0.022231 *  
## pat.karno -1.091e-02  9.891e-01  8.141e-03 -1.340 0.180160    
## meal.cal   2.602e-06  1.000e+00  2.677e-04  0.010 0.992244    
## wt.loss   -1.671e-02  9.834e-01  7.911e-03 -2.112 0.034647 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##           exp(coef) exp(-coef) lower .95 upper .95
## inst         0.9701     1.0308    0.9455    0.9954
## age          1.0129     0.9873    0.9895    1.0369
## sex          0.5674     1.7623    0.3824    0.8420
## ph.ecog      2.4778     0.4036    1.5523    3.9552
## ph.karno     1.0269     0.9738    1.0038    1.0506
## pat.karno    0.9891     1.0110    0.9735    1.0051
## meal.cal     1.0000     1.0000    0.9995    1.0005
## wt.loss      0.9834     1.0169    0.9683    0.9988
## 
## Concordance= 0.648  (se = 0.03 )
## Likelihood ratio test= 33.7  on 8 df,   p=5e-05
## Wald test            = 31.72  on 8 df,   p=1e-04
## Score (logrank) test = 32.51  on 8 df,   p=8e-05

Fitting the model and plotting

cox <- survfit(cox_mod)

plot(cox, main = "Cox proportional hazard model", xlab = "Number of days", 
     ylab = "Prob of survival")

More advanced example of survival models

Source: https://rviews.rstudio.com/2017/09/25/survival-analysis-with-r/

Loading the data and packages

library(survival)
library(ranger)
library(ggplot2)
library(dplyr)
library(ggfortify) ## for better survival plots

data(veteran)
## Warning in data(veteran): data set 'veteran' not found
head(veteran)
##   trt celltype time status karno diagtime age prior
## 1   1 squamous   72      1    60        7  69     0
## 2   1 squamous  411      1    70        5  64    10
## 3   1 squamous  228      1    60        3  38     0
## 4   1 squamous  126      1    60        9  63    10
## 5   1 squamous  118      1    70       11  65    10
## 6   1 squamous   10      1    20        5  49     0

Kaplan Meier Analysis

km <- with(veteran, Surv(time, status))
head(km, 80)
##  [1]  72  411  228  126  118   10   82  110  314  100+  42    8  144   25+  11 
## [16]  30  384    4   54   13  123+  97+ 153   59  117   16  151   22   56   21 
## [31]  18  139   20   31   52  287   18   51  122   27   54    7   63  392   10 
## [46]   8   92   35  117  132   12  162    3   95  177  162  216  553  278   12 
## [61] 260  200  156  182+ 143  105  103  250  100  999  112   87+ 231+ 242  991 
## [76] 111    1  587  389   33

Fitting model with estimates at 1, 30, 60, and 90 days

km_fit <- survfit(Surv(time, status) ~ 1, data = veteran)
summary(km_fit, times = c(1, 30, 60, 90*(1:10)))
## Call: survfit(formula = Surv(time, status) ~ 1, data = veteran)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1    137       2    0.985  0.0102      0.96552       1.0000
##    30     97      39    0.700  0.0392      0.62774       0.7816
##    60     73      22    0.538  0.0427      0.46070       0.6288
##    90     62      10    0.464  0.0428      0.38731       0.5560
##   180     27      30    0.222  0.0369      0.16066       0.3079
##   270     16       9    0.144  0.0319      0.09338       0.2223
##   360     10       6    0.090  0.0265      0.05061       0.1602
##   450      5       5    0.045  0.0194      0.01931       0.1049
##   540      4       1    0.036  0.0175      0.01389       0.0934
##   630      2       2    0.018  0.0126      0.00459       0.0707
##   720      2       0    0.018  0.0126      0.00459       0.0707
##   810      2       0    0.018  0.0126      0.00459       0.0707
##   900      2       0    0.018  0.0126      0.00459       0.0707

Plotting

autoplot(km_fit)

Survival curves by treatment

km_trt_fit <- survfit(Surv(time, status) ~ trt, data = veteran)
autoplot(km_trt_fit)

Plotting patients by age

vet <- mutate(veteran, AG = ifelse((age < 60), "LT60", "OV60"),
              AG = factor(AG),
              trt = factor(trt, labels = c("standard", "test")),
              prior = factor(prior, labels = c("NO", "Yes")))

km_AG_fit <- survfit(Surv(time, status) ~ AG, data = vet)
autoplot(km_AG_fit)

Cox analysis

Fitting model

cox <- coxph(Surv(time, status) ~ trt + celltype + karno + diagtime + age
             + prior, data = vet)
summary(cox)
## Call:
## coxph(formula = Surv(time, status) ~ trt + celltype + karno + 
##     diagtime + age + prior, data = vet)
## 
##   n= 137, number of events= 128 
## 
##                         coef  exp(coef)   se(coef)      z Pr(>|z|)    
## trttest            2.946e-01  1.343e+00  2.075e-01  1.419  0.15577    
## celltypesmallcell  8.616e-01  2.367e+00  2.753e-01  3.130  0.00175 ** 
## celltypeadeno      1.196e+00  3.307e+00  3.009e-01  3.975 7.05e-05 ***
## celltypelarge      4.013e-01  1.494e+00  2.827e-01  1.420  0.15574    
## karno             -3.282e-02  9.677e-01  5.508e-03 -5.958 2.55e-09 ***
## diagtime           8.132e-05  1.000e+00  9.136e-03  0.009  0.99290    
## age               -8.706e-03  9.913e-01  9.300e-03 -0.936  0.34920    
## priorYes           7.159e-02  1.074e+00  2.323e-01  0.308  0.75794    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                   exp(coef) exp(-coef) lower .95 upper .95
## trttest              1.3426     0.7448    0.8939    2.0166
## celltypesmallcell    2.3669     0.4225    1.3799    4.0597
## celltypeadeno        3.3071     0.3024    1.8336    5.9647
## celltypelarge        1.4938     0.6695    0.8583    2.5996
## karno                0.9677     1.0334    0.9573    0.9782
## diagtime             1.0001     0.9999    0.9823    1.0182
## age                  0.9913     1.0087    0.9734    1.0096
## priorYes             1.0742     0.9309    0.6813    1.6937
## 
## Concordance= 0.736  (se = 0.021 )
## Likelihood ratio test= 62.1  on 8 df,   p=2e-10
## Wald test            = 62.37  on 8 df,   p=2e-10
## Score (logrank) test = 66.74  on 8 df,   p=2e-11

Plotting

cox_fit <- survfit(cox)
autoplot(cox_fit)

Notes on cox regression models

  • Need to be cautious when interpreting results because it assumes covariates are constant, and in many cases there are not

  • Hazard: the rate at which events happen. CPH assumes hazards are constant over time

  • Hazard ratio: ratio of hazard between covariates (such as male/female)

  • Look more into Concordance stats if wanting to use ROC curves to assess model performance

Aalen’s additive regression model shows how covariates change over time

aa_fit <- aareg(Surv(time, status) ~ trt + celltype +
                  karno + diagtime + age + prior, data = vet)

aa_fit
## Call:
## aareg(formula = Surv(time, status) ~ trt + celltype + karno + 
##     diagtime + age + prior, data = vet)
## 
##   n= 137 
##     75 out of 97 unique event times used
## 
##                       slope      coef se(coef)      z        p
## Intercept          0.083400  3.81e-02 1.09e-02  3.490 4.79e-04
## trttest            0.006730  2.49e-03 2.58e-03  0.967 3.34e-01
## celltypesmallcell  0.015000  7.30e-03 3.38e-03  2.160 3.09e-02
## celltypeadeno      0.018400  1.03e-02 4.20e-03  2.450 1.42e-02
## celltypelarge     -0.001090 -6.21e-04 2.71e-03 -0.229 8.19e-01
## karno             -0.001180 -4.37e-04 8.77e-05 -4.980 6.28e-07
## diagtime          -0.000243 -4.92e-05 1.64e-04 -0.300 7.65e-01
## age               -0.000246 -6.27e-05 1.28e-04 -0.491 6.23e-01
## priorYes           0.003300  1.54e-03 2.86e-03  0.539 5.90e-01
## 
## Chisq=41.62 on 8 df, p=1.6e-06; test weights=aalen
autoplot(aa_fit)

Tuning cox proportional hazards regression

Source: https://www.youtube.com/watch?v=TrS2M5imOt8&t=9s

additional research on JMbayes2 model

Source: https://drizopoulos.github.io/JMbayes2/articles/JMbayes2.html Primary function is jm()

Basic use example

Goal: to assess the strength of the association between the risk of death and levels of serum bilirubin

library(JMbayes2)
pbc2.id$status2 <- as.numeric(pbc2.id$status != 'alive')
CoxFit <- coxph(Surv(years, status2) ~ sex, data = pbc2.id)

Describe patient profiles for biomarker using a linear mixed model

fml <- lme(log(serBilir) ~ year * sex, data = pbc2, random = ~ year | id) 

Linking survival and longitudinal models

jointFit1 <- jm(CoxFit, fml, time_var = "year")
summary(jointFit1)
## 
## Call:
## jm(Surv_object = CoxFit, Mixed_objects = fml, time_var = "year")
## 
## Data Descriptives:
## Number of Groups: 312        Number of events: 169 (54.2%)
## Number of Observations:
##   log(serBilir): 1945
## 
##                  DIC     WAIC      LPML
## marginal    4363.693 5356.479 -3358.346
## conditional 3538.695 3357.098 -1901.875
## 
## Random-effects covariance matrix:
##                     
##        StdDev   Corr
## (Intr) 1.0032 (Intr)
## year   0.1839 0.4006
## 
## Survival Outcome:
##                         Mean  StDev    2.5%  97.5%      P   Rhat
## sexfemale            -0.1529 0.2694 -0.6499 0.3888 0.5636 1.0007
## value(log(serBilir))  1.2507 0.0866  1.0810 1.4238 0.0000 1.0117
## 
## Longitudinal Outcome: log(serBilir) (family = gaussian, link = identity)
##                   Mean  StDev    2.5%   97.5%      P   Rhat
## (Intercept)     0.7238 0.1720  0.3852  1.0617 0.0000 0.9998
## year            0.2673 0.0382  0.1930  0.3445 0.0000 1.0037
## sexfemale      -0.2639 0.1824 -0.6187  0.0887 0.1527 0.9998
## year:sexfemale -0.0887 0.0407 -0.1687 -0.0094 0.0276 1.0046
## sigma           0.3465 0.0065  0.3342  0.3597 0.0000 1.0105
## 
## MCMC summary:
## chains: 3 
## iterations per chain: 3500 
## burn-in per chain: 500 
## thinning: 1 
## time: 7 sec
  • Assume the instantaneous risk of an event at a specific time \(t\) is associated with the value of the linear predictor of the longitudinal outcome at the same point \(t\)

Traceplots:

  • MCMC stands for Markov chain Monte Carlo

  • Traceplots are used to determine how many iteration it will take to summarize the posterior distribution (convergence) (think Bayes and conditional probabilities)

    • If model converges, then traceplot moves around the mode of distribution
ggtraceplot(jointFit1, "alphas")

Density plot

ggdensityplot(jointFit1, "alphas")